function [dof_map_weight, dof_map_weight_inv] = build_c1_weight2(V, T, E, TE, d)
% function [dof_map_weight, dof_map_weight_inv] = build_c1_weight2(V, T, E, TE, d)
% 
% this version simplify the edge d.o.f
% To calculate the d.o.f. weight for the Argyris elements with degree d
%

% calculate the barycentric coordinate of x- and y- axis 
% with respect to all triangles
[taux1,taux2,taux3] = directional_bary(V,T,1,0);
[tauy1,tauy2,tauy3] = directional_bary(V,T,0,1);

% calculate the edges' normal direction's barycentric coordinate
nE = size(E,1);
normal = zeros(2,nE);
for eg = 1:nE
    dir = (V(E(eg,2),:) - V(E(eg,1),:))';
%    normal(:,eg) = [ 0 -1 ; 1 0]*(dir/norm(dir,2));
    normal(:,eg) = [ 0 -1 ; 1 0]*dir;
end

% prepare space for dof map weight
m = (d+1)*(d+2)/2;
ndofv = 6; 
ndofe = 2*d - 9;
ndoft = (d-4)*(d-5)/2;
n = 3*ndofv + 3*ndofe + ndoft;

nT = size(T,1);
dof_map_weight = zeros(n,m,nT); 
dof_map_weight_inv = zeros(n,m,nT);

%% 6 d.o.f. per vertex
% vertex 1
    
  dof_map_weight(1,1,:) = 1;
  
  dof_map_weight(2,1,:) = taux1; dof_map_weight(3,1,:) = tauy1;
  dof_map_weight(2,2,:) = taux2; dof_map_weight(3,2,:) = tauy2;
  dof_map_weight(2,3,:) = taux3; dof_map_weight(3,3,:) = tauy3;
  
  dof_map_weight(4,1,:) = taux1.*taux1;
  dof_map_weight(4,2,:) = 2*taux1.*taux2;
  dof_map_weight(4,3,:) = 2*taux1.*taux3;
  dof_map_weight(4,4,:) = taux2.*taux2;
  dof_map_weight(4,5,:) = 2*taux2.*taux3;
  dof_map_weight(4,6,:) = taux3.*taux3;

  dof_map_weight(5,1,:) = taux1.*tauy1;
  dof_map_weight(5,2,:) = taux1.*tauy2 + tauy1.*taux2;
  dof_map_weight(5,3,:) = taux1.*tauy3 + tauy1.*taux3;
  dof_map_weight(5,4,:) = taux2.*tauy2;
  dof_map_weight(5,5,:) = taux2.*tauy3 + taux3.*tauy2;
  dof_map_weight(5,6,:) = taux3.*tauy3;
  
  dof_map_weight(6,1,:) = tauy1.*tauy1;
  dof_map_weight(6,2,:) = 2*tauy1.*tauy2;
  dof_map_weight(6,3,:) = 2*tauy1.*tauy3;
  dof_map_weight(6,4,:) = tauy2.*tauy2;
  dof_map_weight(6,5,:) = 2*tauy2.*tauy3;
  dof_map_weight(6,6,:) = tauy3.*tauy3;
  
% vertex 2

tau = taux1; taux1 = taux2; taux2 = taux3; taux3 = tau;
tau = tauy1; tauy1 = tauy2; tauy2 = tauy3; tauy3 = tau;

  dof_map_weight(7,16,:) = 1;
  
  dof_map_weight(8,16,:) = taux1; dof_map_weight(9,16,:) = tauy1;
  dof_map_weight(8,17,:) = taux2; dof_map_weight(9,17,:) = tauy2;
  dof_map_weight(8,11,:) = taux3; dof_map_weight(9,11,:) = tauy3;
  
  dof_map_weight(10,16,:) = taux1.*taux1;
  dof_map_weight(10,17,:) = 2*taux1.*taux2;
  dof_map_weight(10,11,:) = 2*taux1.*taux3;
  dof_map_weight(10,18,:) = taux2.*taux2;
  dof_map_weight(10,12,:) = 2*taux2.*taux3;
  dof_map_weight(10,7,:) = taux3.*taux3;

  dof_map_weight(11,16,:) = taux1.*tauy1;
  dof_map_weight(11,17,:) = taux1.*tauy2 + tauy1.*taux2;
  dof_map_weight(11,11,:) = taux1.*tauy3 + tauy1.*taux3;
  dof_map_weight(11,18,:) = taux2.*tauy2;
  dof_map_weight(11,12,:) = taux2.*tauy3 + taux3.*tauy2;
  dof_map_weight(11,7,:) = taux3.*tauy3;
  
  dof_map_weight(12,16,:) = tauy1.*tauy1;
  dof_map_weight(12,17,:) = 2*tauy1.*tauy2;
  dof_map_weight(12,11,:) = 2*tauy1.*tauy3;
  dof_map_weight(12,18,:) = tauy2.*tauy2;
  dof_map_weight(12,12,:) = 2*tauy2.*tauy3;
  dof_map_weight(12,7,:) = tauy3.*tauy3;

% vertex 3

tau = taux1; taux1 = taux2; taux2 = taux3; taux3 = tau;
tau = tauy1; tauy1 = tauy2; tauy2 = tauy3; tauy3 = tau;

  dof_map_weight(13,21,:) = 1;
  
  dof_map_weight(14,21,:) = taux1; dof_map_weight(15,21,:) = tauy1;
  dof_map_weight(14,15,:) = taux2; dof_map_weight(15,15,:) = tauy2;
  dof_map_weight(14,20,:) = taux3; dof_map_weight(15,20,:) = tauy3;
  
  dof_map_weight(16,21,:) = taux1.*taux1;
  dof_map_weight(16,15,:) = 2*taux1.*taux2;
  dof_map_weight(16,20,:) = 2*taux1.*taux3;
  dof_map_weight(16,10,:) = taux2.*taux2;
  dof_map_weight(16,14,:) = 2*taux2.*taux3;
  dof_map_weight(16,19,:) = taux3.*taux3;

  dof_map_weight(17,21,:) = taux1.*tauy1;
  dof_map_weight(17,15,:) = taux1.*tauy2 + tauy1.*taux2;
  dof_map_weight(17,20,:) = taux1.*tauy3 + tauy1.*taux3;
  dof_map_weight(17,10,:) = taux2.*tauy2;
  dof_map_weight(17,14,:) = taux2.*tauy3 + taux3.*tauy2;
  dof_map_weight(17,19,:) = taux3.*tauy3;
  
  dof_map_weight(18,21,:) = tauy1.*tauy1;
  dof_map_weight(18,15,:) = 2*tauy1.*tauy2;
  dof_map_weight(18,20,:) = 2*tauy1.*tauy3;
  dof_map_weight(18,10,:) = tauy2.*tauy2;
  dof_map_weight(18,14,:) = 2*tauy2.*tauy3;
  dof_map_weight(18,19,:) = tauy3.*tauy3;

for k = 1:size(T,1)
%% 2*d - 9 d.o.f. per edge
% edge 1: v2-v3
    eg = TE(k,1);
    [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
    dof_map_weight(19,[13 18 19],k) = [tau1 tau2 tau3]; 
% edge 2: v3-v1
    eg = TE(k,2);
    [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
    dof_map_weight(20,[9 10 6],k) = [tau1 tau2 tau3]; 
% edge 3: v1-v2;
    eg = TE(k,3);
    [tau1,tau2,tau3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), normal(1,eg), normal(2,eg));
    dof_map_weight(21,[8 4 7],k) = [tau1 tau2 tau3];  

%% There is no identity matrix for the rest inner d.o.f.

    idx = [1 2 3 4 5 6 16 17 11 18 12 7 21 15 20 10 14 19 8 9 13];
    dof_map_weight_inv(idx,:,k) = inv(dof_map_weight(:,idx,k));

end

end

function [tau1,tau2,tau3] = bary(V1,V2,V3,x,y)
x1 = V1(1); x2 = V2(1); x3 = V3(1);
y1 = V1(2); y2 = V2(2); y3 = V3(2);

J  = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
J1 = (x2 - x )*(y3 - y ) - (x3 - x )*(y2 - y );
J2 = (x  - x1)*(y3 - y1) - (x3 - x1)*(y  - y1);
J3 = (x2 - x1)*(y  - y1) - (x  - x1)*(y2 - y1);

% barycetric for (0,0);
J10 = x2.*y3 - x3.*y2;
J20 = (x3 - x1).*y1 - (y3 - y1).*x1;
J30 = (y2 - y1).*x1 - (x2 - x1).*y1;

tau1 = (J1 - J10)./J;
tau2 = (J2 - J20)./J;
tau3 = (J3 - J30)./J;
end


